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1. INTRODUCTION 



In 1981-1983 a computer program was developed which solves numerically 
two-dimensional Euler equations using the Godonov method (Ref. 1). The 
program was intended primarily for turbomachinery application, but with 
modification to the boundary conditions, application to any problem which is 
modeled by the Euler equations is possible. 

This report covers one year of extensive testing of the code for 
different flow regimes. This testing was done with three main objectives: 

a) To gain experience in simulation of the various problems related to 
the research activity in the Turbopropulsion Laboratory, 

Naval Postgraduate School, Monterey, California. 

b) To test code accuracy by comparing numerical with analytical 
solutions . 

c) To find the source of errors for various nur\erical simulations and the 
basic Godunov method in order to reduce these errors. 

Some of the results of these experiments with the code were reported in 
other papers and appear as appendices to this report. Other cases were 
simulated, but not published, and will be reported here. Some unsuccessful 
attempts to improve the code accuracy are also described. 

The topics studied and the progress made are described in the following 
sections . 

1) The problem of the gradual opening of wave rotor passage. 

2) Numerical modeling of the nonsteady thrust produced by intermittent 
pressure rise in a diverging channel. 

3) Numerical techniques for wave rotor cycle analysis (A. Mathur). 

4) Development of a grid generation program in order to generate a grid 
over wedge-arc configuration. 

Test case of the supersonic flow in a channel with arc, wedge, 
wedge-arc, sinsoidal arc-bumps. 

6) Modification of the Godunov method to include shock fitting using the 
oblique shock wave theory. 

7) Development of the one-dimensional code based on Verhoff’s method. 

8) Modification of the program for cylindrical geometry problems, .and 
modeling of the flow in CDTD. 

9) ModeLing of the flow in a detonation engine. 

10) Development of the SELCO Method. 

11) Grid generation for CDTD inlet 
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2 . 



PROBLEMS CONSIDERED 



2 . 1 The Gradual Opening in Wave Rotor Passage . 



Influence on the flow pattern of the gradual passage opening in the wave 
rotor was studied on the basis of a numerical simulation. It was found that, 
in most cases, a significant volume of the passage will have rotational flow, 
which should lead to the mixing between the driver and driven gases. In some 
cases, losses will occur as a result of the multiple reflections of the shock 
and pressure waves from the passage walls. . It was shown that the interface 
between driven and driver gas will be oblique to the passage walls when the 
passage opens gradually, and that the interface can -retain its obliqueness to 
the walls. The results for the rectangular and skewed passages are reported 
in Appendix A. 



In the case of the skewed passage, the perfect compression of the driven 
gas with very small mixing losses could be achieved if the channel opens 
according to a certain formula of the optimum opening velocity. It was found 
that the velocity of opening could be predicted by the equation: 



V 



open 



^sh 

sT n"~ L sk 



( 1 ) 



where: ^open = the velocity of the gradual opening 

v sh = the shock wave velocity in the driven gas 

■^sk = t ^ ie of the skewed passage 



2 . 2 Numerical Modeling of the Nonsteady Thrust Produced by Intermittent 
Pressure Rise in a Diverging Channel. 



The dynamics of the expansion of detonation products In a diverging 
channel were investigated numerically. The influence of particular features 
in the expansion process, such as the presence of reversed flow and 
propagation of hammer shocks on the production of thrust were examined. 
Sequential expansions of detonation products were also modeled and It was 
concluded that in order to maintain a high frequency periodic mode of 
operation for propulsion applications, the channel should be refilled with 
ambient air after each expansion. The influence of the ratio of ambient air 
to detonation products volume on the dynamics of the thrust production and on 
the impulse generated during the expansion are also reported. The 
presentation of the results of this investigation is given in Appendix B. 

2 . 3 Numerical Techniques for Wave Rotor Cycle Analysis. 



With the author’s participation, A. Mathur developed a one-dimensional 
code based on Godunov’s method for wave rotor cycle ina lysis. Because the 
resolution of the code on the interfaces was poor, it was decided to develop a 
one-dimensional code based on the random choice method. Non-standard boundary 



conditions for this problem were developed. A paper was written and accepted 
for publication in the proceedings of the ASME Symposium on Nonsteady Flows. 
The paper is attached as Appendix C. 

2 . 4 Development of the Grid Generation Program for Wedge-Arc Configurations . 

The grid generation routine, which was used for the cascade and wave 
rotor simulations did not allow a full control of the grid parameters. This 
control was especially needed for the simulation of the test problems, since 
we wanted to insure that the grids have the same general parameters for the 
different wedge-arc configurations considered. 

A program which allows the grid generation with full control over grid 
parameters was developed. The program, with an example of the grid, is given 
in Appendix D was delivered to the Naval Postgraduate School and stored in an 
archive file. 

2 . 5 Test Case of the Supersonic Flow in a Channel with Arc, hedge, Wedge-Arc, 
Sinusoidal Arc-Bumps . 



Performance of the Godunov code was tested for numerous supersonic and 
subsonic internal flow problems. The tasks included: 

a) An attempt to determine what causes the nonsyme trical behavior for 
subsonic flow case and whether it is possible to increase the 
accuracy of the basic Godunov code by changing boundary conditions. 

b) A comparison of the numerical solution with analytical predictions 
for supersonic flow. 

For subsonic flows it was determined that the main source of errors was slop 
discontinuity at the corners of the "bump" . When the circular arc bump was 
replaced by a sinusoidal arc “bump" , the symmetry of the subsonic solution 
improved substantially. This resulted from the fact that in the case of the 
sinusoidal arc only the second derivatives are discontinuous at the 
surface. It was found that the boundary conditions at the surface, which are 
used by the basic Godunov method (solution of the Riemann problem between the 
physical point and it mirror image with the same pressure and density but with 
the negative velocity), produce the most consistent results for various 
problems. Results for this simulation are shown in Figure 1. The grid is 
shown in Figure 2. 

Efforts to improve the subsonic flow solution were discontinued when it 
was realized that a more objective criteria than symmetry of solution was 
needed. For this reason, work on supersonic flows was begun, since, for these 
conditions, some analytical solutions are available. First, the method for 
supersonic flow (M = 2.2) in a channel with a 10% thick circular arc "bump was 
tried. The grid for this case is shown in Figure 3. In Figure 4 we show 
results for this geometry when the inlet Mach number is 2.0. Results are 
shown in the form of the distributions of the pressure coefficient over the 
surface of the lower wall of the channel. Numerical results correspond very 
well to the analytical In this case. It was found, however, that when other 
flow parameters are compared with the analytical solution, there was 
substantial disagreement between the resits. Figure 5 shows the comparison of 
the pressure coefficient, entropy, density, velocity and Mach number that were 



3 



calculated numerically with the analytical solution. The comparison is done 
for supersonic flow with '1=1.6 over 5% thick wedge-arc "bump". As shown in 
Figure 5, the errors are produced mostly on the shock waves and that the 
largest is the relative error in entropy. The fact that only pressure is 
simulated accurately was confirmed for the numerous test cases in which the 
arc thickness and the wedge angle varied in the wide range of values. 

2 . 6 Modification of the Godunov Method based on the Shock Fitting using the 
Oblique Shock Wave Theory . 

A shock fitting was attempted to reduce the errors on the oblique shocks. 
The procedure for fitting was as follows: 

1) Calculate the values using Godunov's method. 

2) Using oblique shock wave theory, calculate the exact value of the 
parameters. This calculation is based on the assumption that the 
pressure is computed correctly by the Godunov method. 

3) Substituting values behind the shock for the exact values. It was 
found that the procedure did not work if it was applied as 
described above. Many modifications of this prodedure were tried. 
Some involved interpolation of the values ahead and behind the 
shock. The last approach did improve the accuracy of the entropy 
calculation but introduced oscillations. 

2 . 7 Development of the One-Dimensional Code based on the Verhoff Method . 

A new formulation of the Euler equations and numerical method to solve 
them were proposed by Verhoff (Ref. 2). A basic code implementing his method 
for one -dimensional flow was developed. In Figure 6, results are shown for 
the case of two colliding shock waves in a tube. These results are a 
reconstruction of the results presented in Ref. 1 in Figure 5. Comparison of 
these two figures shows that our code implements the basic Verhoff method 
correctly. The listing of the program is presented in Appendix E. 

We were unable to implement the shock fitting in our code. Comparison 
between the analytical solution and numerical was 10-15% in error for 
velocity, pressure and density in the region of the shock wave, if the shock 
was not fitted. 

2.8 Modification of the Program for Cylindrical Geometry Problems, and 
Modeling of the Flows in CDTD. " 



The Godunov code was modified for simulation of the flow with 
cylindrical velocity. The modified program was used for modeLing of the flow 
field in CDTD device (Ref. 3). The task of the modeling was to test how the 
rotational component of velocity will influence flow field. 

In Figure 7, one of the geometries of the CDTD is shown with a 
computational grid covering the domain of integration. In Figure 8, the 
Iso-Mach lines for the geometry are shown for the inlet flow conditions 
indicated in figure 7. 



The conclusions of preliminary research was that the program models the 
flow correctly and that, for the tested conditions , the results are close Co 
those given by the program developed by Hirsch as discussed in Ref. 4. 

2 . 9 Modeling of the Flow in Detonation Engine . 

The flow field in the projected combustion chamber of a detonation 
engine was simulated numerically using the Godunov code. The cask of chis 
research was to study the dynamics of the discharge of the high pressure gas 
after detonation. A number of configurations were tried. In Figure 9, an 
example of this type of simulation can be seen. Pressure in the detonation 
chamber drops from 15 atm to 0.55 atm in 0.00131 second. It was important to 
establish that time since it is one of the limitations on the deconation 
frequency. In this example, the ability of the program Co simulate large 
gradients of pressure can be seen. 

2.10 Development of the SELCO Method . 

A new method for supersonic flow simulation was developed. A full 
description of the method, with illustrations of its application to the 
supersonic flow over the wedge, is given in Appendix D. This method 
demonstrates that inaccurate modeling of the oblique shock waves produced by 
the Godunov method is the result of the obliqueness of the shock wave with 
respect to the edges of the cells of the computational grid covering domain of 
integration. It was shown that only when che shock surface is parallel to the 
two opposite edges of the cell that the oblique shock can be accurately 
calculated . 

A new method of the Local Cell Orientation, - SELCO, was proposed to 
allow local reorientation of the cells in the vicinity of the shock waves. 

The efficiency of the SELCO method was demonstrated for the simulation of the 
oblique shock waves in the supersonic flow over the wedge. 

2.11 Grid Generation for CDTD Inlet. 



For a given distribution of the pressure measurement ports on the CDTD 
blade surface, a computational grid was developed. Several variants of the 
CDTD geometry were implemented. However all grids which were generated proved 
to lead to large computational errors. This work was not continued due to time 
constraints on the project. 

3. Conclusions 



The principal findings of the program of code cescing were as follows: 

1. The code has a robust ability to analyse unsteady time dependent 
flows and provide qualitative information on complex shock patterns 
and pressure fields. 

2. With attention to boundary conditions and particularly grid 
structure it is possible to reproduce analytical flow solutions with 
great accuracy. 
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3. The solutions produced are sensitive to grid orientation. The grid 
must be re-orientated locally along shock or discontinuity planes. 



4. Care must be used in selecting the grid distribution on areas of 
rapid curvature change. 

Despite the program sensitivity under certain conditions, with skilled 
use it can be most instructive in establishing wave and shock structures 
in complex flow situations. A significant example is the unsteady 
boundary of Section 2.1 where in it was possible to generate a minimum 
loss passage angle for a wave rotor device. 
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Figure 1. Pressure coefficient and surface Mach 

number for the subsonic flow (M . = 0.5) 
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Figure 3. Grid for circular arc bump. 
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Mach number for the flow in the channel at M 
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15 



APPENDIX A 



The Problem of Gradual Opening in Wave Rotor Passages Gradual Opening of 
Skewed Passages in Wave Rotors • 



16 



The Problem of Gradual Opening in Wave Rotor Passages 



S. Eidelman* 

Naval Postgraduate School , Monterey, California 



The influence on the flow pattern of (he gradual passage opening in the wave rotor is studied on the hasis of 
numerical simulation. It is found that, in most cases, significant volume of the passage will hu\e rotational flow, 
which should lead to the mixing between the driver and driven gases. In some cases, losses will occur also as a 
result of the multiple reflection of the shock and pressure waves from the passage walls. Il is >hown (hat the 
interface helwcen driven and driver gas will he ohliijue to the passage walls, when the passage opens gradually 
and (lie interface can retain its obliqueness to the walls. 



Introduction 

W AVE rotors are devices that use wave propagation 
through the fluid in a rotor passage to transfer energy 
from one fluid to another 1 or to transfer energy from one 
fluid to rotor shaft and another fluid. 2 Wave rotors, wave 
engines, wave pressure exchangers, wave equalizers, and 
Cornprex 1 are all based on the same idea of energy exchange 
in the unsteady waves, and the topic of the present study is 
relevant to all of these devices. 

Wave rotors offer the potential for significant im- 
provements in air-breathing engine propulsion cycles because 
they can be self-cooling, since both high-pressure gas and cold 
low-pressure gas use the same passage for alternate periods of 
time in the cycle. Combining a wave rotor with conventional 
turbomachinery components shows promise of significant 
reduction in specific fuel consumption without weight or size 
penalties. 

The principles of operation of wave rotor devices and their 
commerical applications can be found in Refs. 1-3. For 
completeness, the general scheme of a wave pressure ex- 
changer will be described, as illustrated in Fig. 1. One gas 
(driver) at high pressure is used to compress a second gas 
(dmen). The process is arranged to occur in tube-like 
passages with trapezoidal cross section located on the per- 
iphery of a rotating drum or rotor. The compression is 
achieved successively in each rotor passage by means of 
compression waves or shock waves generated by the entering 
driver gas. The compressed driven gas is drawn off from the 
end of the passage when it aligns with an outlet port. The 
driver gas then undergoes a series of expansions to a lower 
pressure and is scavenged out by freshly inducted driven gas at 
approximately the same pressure level. This fresh “charge” is 
then compressed by the high-pressure driver gas and the cycle 
repeats itself. Steady rotation of the drum sequences the ends 
of the passages past stationary inlet ports, outlet ports, and 
cndwalls. This establishes unsteady but periodic flow pro- 
cesses within the rotating passages and essentially steady flow 
in the inlet and outlet ports. The passage may be oriented 
axially as in Tig. I or at a stagger angle. In general, wave 
machines used as pure pressure exchangers (e.g., “Coin- 
prex” 1 ) have usually axially oriented passages, while those 
with staggered passages may drive a shall, since shaft work 
extraction is possible with this latter configuration. 
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In the design of the wave rotor, it is very desirable u> 
determine the optimal ratio between the width and length ot 
the single passage of the rotor. Usually only two parameters 
are evaluated to determine this ratio; skin friction losses and 
bypass losses. The number of passages on the rotor should be 
minimal to minimize skin friction losses. On the other hand, 
the passage should be narrow compared to the port widths to 
reduce flow bypass between inlet or outlet ports. The transient 
process of the passage opening or closing (as the passage end 
moves across a port or moves from a port to be closed b> the 
endwall, respectively) usually is not considered in the design. 
It is generally assumed that the passage opening or closure 
occurs instantaneously. 

Pearson 2 tried to take into account the gradual opening ot 
the passage, assuming that the air in the passage is com pres >ed 
in a series (usually three) of discrete compression waves which 
converge and ultimately merge to form a shock wave. This 
allowed him to design a complicated wave machine c>cle im a 
rotor using relatively short passages. However, since t lie 
technique was one-dimensional it could not reveal t lie 
peculiarities of this essentially two-dimensional flow and 
would be valid only for very weak waves. 

In the present study, by means of numerical modeling, we 
will examine in real-time how the gradual opening mlluences 
the wave formation in the wave rotor passages and how that 
should affect the rotor design. 

Model 

The assumptions involved in the numerical simulations ate 
described as follows. 

The flow in each passage of the wave rotor is unsleadv and 
periodic. At the same time, the flow through t he port-, is 
(ideally) steady. 1 2 The peripheral width ot the port is usu.dlv 
equal to several widths ot a single passage, and herein it is 
assumed that the How in the inlet or outlet port remains 
stationary when the passage-end encounters the port. For ihis 
reason the region of the port is not included in the miu- 
piitai ionat domain shown in Tig. 2. 

The tiine-dependem process of the passage-end translaimg 
across the region of the inlet or outlet port will be referred as 
the “gradual opening” of the passage. The passage opening 
process will he referred as “instantaneous.” when the as- 
sumption is made that the passage instantaneously opens io 
the port area and is subjected immediately to the stead) I low 
conditions at the port. 

It is assumed Mini the flow is in viscid and can be modeled by 
the linler equations. Ihe unsteady two-dimensional Euler 
equations cun he written in conservation law form as 

du op do 

7 + d\ + 7Y =0 
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where 
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pu 




pr 




pu 




p + pu 2 




plW 


u= 


pv 


F= 


pill’ 


G = 


P + pv : 




^ J 




, (e + p) ii , 




, (e + p)v , 



(I) 

Here p is the density, u and v the velocity components in the X 
and Y coordinate directions, p the pressure, and 7 the ratio of 
specific heats. T he energy per unit of volume, e, is defined by 



where e = P( ( 7 - I)p is the internal energy. We look for the 
solution of the system of equations represented by Eq. ( 1 ) in 
the computational domain shown in Fig. 2 in time 1 with the 
following conditions at the domain boundaries: a) solid wall 
along segments 1-3 and 2-4, b) outflow along segment 3-4, 
and c) inlet along segment 1 - 2 . 

It is assumed that initially at time / = 0 , the passage of the 
wave rotor is filled with stationary gas at ambient conditions. 
When instantaneous opening of the wave rotor passage was 
simulated, it was assumed that at time / = 0 , the flow at the 
inlet 1-2 was equal to the steady flow in the port. When 
gradual opening of the passage was simulated, at time / = 0 , 
the inlet was closed and solid wall boundary conditions were 
imposed at the inlet 1-2. Then, this bo&ndary condition was 
gradually replaced by the flow condition at the inlet port. The 
length uncovered to the inlet port region, where the solid wall 
boundary conditions were replaced by the inlet port con- 
ditions, was determined using the elapsed time and the 
velocity of the passage relative to the inlet port. 

The Godunov method was used to obtain a numerical 
solution of Eq. ( 1 ) with the described boundary and initial 
conditions. Details of the implementation of the method and 
boundary conditions are given in Ref. 4. 

The flowfield was modeled for the rectangular passage with 
a width of 0.02 m and a length of 0.12 m. The grid covering 
the computational domain of the passage is shown in Fig. 2. 

Results and Discussion 

The following initial conditions were assumed for the air 
initially in the passage: 

P () -] atm p 0 - 1.2 kg/m 3 U 0 = 0 V () = 0 

The driver gas entering through the port at the left-hand end 
was assumed to have the following properties: 

/T, = l. 8 atm p (/ = 1 .81 kg/m 3 U d = 1 50 m/s V d - 0 

These conditions correspond to a practical situation in a wave 
rotor when a passage filled with a quiescent fresh charge of air 
at ambient conditions encounters an inlet port supplying hot, 
high-pressure driver gas. 

If the assumption is made that the passage instantaneously 
opens and is subjected to the conditions of the inlet flow at the 
left boundary 1-2 (see Fig. 2), then a perfectly one- 
dimensional flow pattern should develop in the passage. The 
results of the modeling of these conditions are presented in the 
form of pressure contours at progressively larger time steps in 
Fig. 3. Time /-() corresponds to the moment when the 
passage opens. Instantaneous opening of the passage is seen in 
Fig. 3 to lead to an immediate formation of a shock wave and 
subsequent propagation in the passage from the left to the 
right. The flow conditions at the inlet port arc matched to the 
parameters of the rarefaction wave which effectively will 
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oils passage opening. 



cancel the rarefaction wave moving toward the inlet at the 
end of the passage. For this reason the flow in the passage 
no additional discontinuity. This situation is typical to 
wave rotor, where flow conditions at the ports are ehosei 
such a way that waves do not propagate from the passe 
into inlet or outlet ports. The llow in the ports will (heret 
remain steady. In the ease shown in lig. 3, the shock \\ 
was found to propagate with a velocity, = 446 m/s, 
interface with a velocity f = 150 m^s, and all the parame 
examined were confirmed accurately using onc-dimensu 
gasdynamics relationships. 

Most of the approaches used in wave rotor design are ba 
on the assumption that the waves are one-dimensiona 
nature. When the gas is compressed by a weak shock w; 
for instance, assumptions are made that, I) the proees 
isoentropic, 2 ) the hot and cold gases are strictly separatee 
a planar interface, and 3) the flow is everywhere irrotatio 
This leads to the very high efficiencies projected tor 
compression. If the passage is wide enough so that \ise 
effects can be neglected, this model of the compression in 
wave rotor passage is tealistic, but only tor insiaiitane 
passage opening or for a very long passage. Results presei 
below illustrate how the gradual passage opening al feels 
compression process. 
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In Fig. 4, the pressure and Much number contours in the 
passage are shown at a sequence of times lor the case ot a 
gradual opening to t he inlet port with a velocity of 100 m/s. 
Parameters for the gas in the passage, belore opening begins, 
and in the inlet port are the same as lor the case of in- 
stantaneous opening. The dynamics of the How development 
seen in Fig. 4a is very different from that shown in Fig. 3. 
First, curved pressure waves radiate Iroin the uiitial small 
opening appearing at the lower corner of the inlet on the left 
side of the passage {( = 0.044 ms in Fig. 4). Subsequently, 
thestT waves reflect from the upper wall of the passage and at 
time / = 0.125 ms have formed a row of compression waves 
which have approximately straight fronts normal to the wall 
of the passage. Initially, compression waves of this kind 
occupy a small part of the region with disturbed gas. In the 
region well behind the front of this quasi-one-dimensional 
propagation, compression waves arc curved and are the result 
of the interaction between the waves reflecting back and forth 
between the lower and upper walls of the passage and the new 
waves created by the progressive opening of the port to the 
passage. Since it is possible to see in Fig. 4 for the times 
/ = 0.084 and 0.166 ms, respectively, the flow behind the 
quasi-one-dimensional region is highly rotational and is 
relatively constant pressure in the A' direction. The passage 
became fully opened only at / = 0. 2 ms. At this time com- 
pression waves are propagating along the length of the 
passage and the pressure rises gradually from 1 atm at the 
right end to 1.8 atm at the left end. The front of converging 
compression waves will eventually become a shock wave, but 
this will occur at some later time and outside the com- 
putational domain. 

It was concluded from the case presented in Fig. 4 that a 
passage length-to-width ratio of 6 will lead to very high 
mixing losses and iionuniform and inefficient compression, 
since in this case the region of rotational flow occupies half of 
the passage volume. Additional losses will be produced 
because of the highly rotational flow within each of the two 
gases. 

In Fig. 5, pressure and velocity contours are shown for the 
case of gradual opening of the passage with a velocity of 200 
m/s. Full opening of the passage in this case occurred at 
/ = ().! ms. A curved shock wave is formed at / = 0.044 ms. 
i This shock wave (see Fig. 5a) partially reflects from the upper 
wall of the passage and then, gradually converging with its 
main Iront, forms an almost planar shock front at the time 
/ = (). 252 ms. However, even then t he flow is highly rotational 
behind the shock front, and the region of rotational How 
occupies one-third of the passage volume. In this region the 
gas velocity increases from M- 0.3 at the upper wall to 
M -0.52 at the lower wall. 

Whet the velocity ol t lie passage opening is increased to 5()t) 
, m/s hce lie. 6) I he passage becomes fully open at ( - 0.04 ms. 
Because of the last opening of the passage, the shock wave at 
t he time t =0.04 ms is less curved and only a small fraction of 
I the shock Iront rolled s from the upper wall. This reflected 
part of the shock front converges with the main front at 
/ = 0.172 ms. From that time on, the How pattern in the 
passage is mostly one-dimensional with a small and 
weakening region of rotational How behind the main front. 
Nevertheless, even for this case, with passage length/ passage 
width 3. there will be very high mixing losses, with a large 
T part ol the passage volume subjected to rotational How. 

1 In outer to study how the strength of the shock wave in- 
fluences the flow pattern developing in the passage, an ad- 
ditional case was simulated where the parameters ol thcdiiver 
, las at the inlet port area were increased; namely, to 

!*,, — 2.85 at m U lf - 283 in/s v,t 4 kg/ in' 

\s in the pievious case, the parameter of the driver gas were 
hosen hi such a way that waves do not propagate back into 
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are shown in Fig. 7. The velocity of the passage opening was 
1 =250 m/s. It was concluded from the results presented in 
Fig. 7 that ail increase in the initial shock strength leads to 
stronger ret lections and substantial increases in the How 
rotation. This, in turn, will lead to increased losses because of 
mixing of the diiver and driven gases. The pattern of multiple 
rel lections of the shock wave from the upper and lower walls 
ol t lie passage can he seen in Fig. 7a. At time t = 0.073 ms the 
shock wave rcllcctcd from the upper wall is propagating 
towaid the lower wall. Fart ot the reflected shock wave Iront 
converges with the main shock wave, and part begins to 
rellect from the lower wall {t =0.109 ms), the multiple 
rel lection weakens the rcllcctcd shock, but the ret lection 
between the walls of the passage continues and can be 
lollowcd tor ( 0.21 1 and 0.245 ms. It is clear I hat lor these 

conditions even passages with a length-to-width ratio ol 6 will 
have substantial mixing losses because ol the rotational I low. 

A dit feieni situation is lound for the passage opening to the 
inlet pori ol the dnveii gas. \i this poll, the pressure and 
velocity ol the (driver) gas in the passage should be matched 
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It die passage opens instantaneously, (he llowtield in the 
passage will have only one discontinuity — the interlace be- 
tween the fresh air and exhaust gas entering and leaving the 
passage, respectively. To model this condition tor the gradual 
passage opening, it was assumed that in the passage: I*,, = I 
ami. />„ 0.5 kg/m\ U, t - 150 m/s, and \'„ 0; in the port: 

/'/ t aim, /i (/ ~t.4 kg/ni’, Uj - 150 m/s, and V ti 0 X he 
velocity ot the passage opening was I -200 nPs. 

Results lor tins simulation are presented in big. K It can be 
seen dim, since at first moment most ol the inlet cross section 



is blocked by the wall, a rarefaction wave reflects from 
inlet wall. Hie passage I nil opening occurs at / = 0. 1 ms. I ; 
this time on the pressure deviation in ihe How field wcu! 
and at / = 0.22 ms the pressure is almost completely unit 
and equal to the static pi cssure m both the tmdisim Kd di 
air ami exhaust gas I igme Sb shows that the imcitace 
carry the imprint ol the gradual passage opening a long 
alter the passage becomes fully open. There is no disMp; 
mechanism included in the mathematical model used m 
study to lorce the interlace to become normal to the pas 
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walls. Therefore, the interface will “remember” the dynamics 
of the gradual passage opening in passages with large lengih- 
to-width ratio. 

Conclusions 

It is shown in the present study, on the basis of numerical 
modeling, that the dynamics of the passage opening sig- 
nificantly affects the flow pattern in the passages of wave 
rotor devices. For rectangular axial passages, even when the 
velocity of the passage opening is 500 m/s, a one-dimensional 
flow pattern forms only in passages with a length-to-width 
ratio larger than 3. In the region preceding the formation of 
the one-dimensional flow pattern, the flow is rotational and, 
in some instances, contains shock and pressures waves which 
repetitively reflect from the upper and lower walls of the 
passage. In practice (his would lead to significant mixing 
between the driver gas (e.g., exhaust gas) and driven gas (e.g., 
fresh air) and reduce the efficiency of the engine cycle. With a 
reduction in the velocity of the passage opening, the volume 
of the mixing region increases as a result of the rotational 
flow which develops from the slow opening of the passage. 
With an increase in shock wave strength, the volume of the 
mixing region increases as a result of the rotational flow 
which develops form multiple reflections of the shock and 
pressure waves from the passage walls. 

Numerical modeling of the process of gradual passage 
opening at the inlet port of the driven gas revealed that the 
interface between the gases will move all the way through the 
passage with a frozen pattern of distortion or obliqueness to 
the passage walls. The interface oblicuieness increases when 
the velocity of the passage opening decreases. 



In all, it can be concluded that taking into account 
gradual passage opening is essential tor the wa\e mac 
design, not only for proper timing and wave arrangement 
also because of the losses which will occur due to mixing 
wave reflections. The gap between projected and achn 
efficiencies of wave machines 2,6 - ■’ could be parualh du 
neglect of the effects of the gradual passage opening w lucl 
studied herein. 
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Gradual' Opening' of' Skewed' Passages' in' Wave Rotors 



Introduction 



The problem of gradual opening of rectangular axial passage in wave 
rotors was studied in our article 1 published in the January issue of Journal 
of Propulsion and Power. There we have defined the problem of gradual 
opening in the wave rotor passage and the mathematical model which was used 
to simulate the opening process. In this article we will assume that the 
reader is familiar with the reference 1 f and it could be regarded as an 
extension to the that reference. 

If the wave rotor is used to produce shaft power, it's passages should 
be skewed in one form or another. In this paper we will analyze the gradual 
opening of a skewed passage and will examine the conclusion drawn for the 
rectangular axial passage for this more general geometry of the passage. 

Pesults and' blscussion 



The main conclusion of the study on gradual opening of rectangular 
passage is that in order to minimize the mixing losses caused by rotational 
flow in the passage, the opening velocity should be very high. In the 
limit, instantaneous opening of the wave rotor passage will lead to one 
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dimentional flow pattern in the passaae will have minimal mixina losses. 

When the passaae of the wave rotor is skewed, even instantaneous openina of 
the passaae will ' not lead to development of the one dimensional flow pattern 
with small mixina losses. That statement will be demonstrated in the 
followina example. 

Let's model the openina process for the passaae 0.02 m wide and 0.24 m 
Iona* It has left and riaht hand inlets parallel to the v-axis and upper 
and lower walls of the passaae form 60° anale with positive direction of the 
x-axis. It is assumed that initially air in the passaae is at the followina 
conditions • 

P = 1 atm; P = 1.2 kq/M 3 ; U = 0? V =0 
O O o o 

The driver aas enterina throuah the port at the left hand end was assumed to 
have the followina properties: 

P, = 1.R atm; P, = 1 .PI kcr/M 3 ; n = 75 M/sec; V = 129.9 M/sec 
d d d d 

The conditions for the driver and driven aas are the same as in case of the 
rectanaular passaae which we reported in the previous paraaraph. 

In Fiaures la and 1h results for simulation of the instantaneous 
openina of the skewed passaae are shown in form of pressure and velocity 
contours at the seauences of times. The flow pattern near the inlet in 
Fioure 1 is hiahly rotational which suaaests very hiah mixina losses. That 
is caused partially by reorientation of the shock wave. At the first 
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moments of time after the openinq of the passaqe the shock wave between the 
driven and driver cases is oblique to the lower and upper walls of the 
passaqe* At later times this shock wave turns and becomes normal to the 
upper and lower walls. Thus, incontrast with rectanqular for skewed 
qeometry passaqes there is no obvious condition for minimizinq the mixinq 
losses caused by the inlet openinq. 

Let's find the conditions for openinq of the skewed passaae which will 
lead to minimal mixinq between driver and driven qases . As we have stated 
before, rotational flow in the instantaneously opened skewed passaqe, will 
develop because of the rotation of the shock wave from the parallel position 
to the inlet (initially) to the position where it is normal to the upper and 
lower walls of the passaqe (when the flow is developed). If we could form a 
shock wave which will be all the time normal to the lower and upper passaae 
walls then the rotational flow and mixinq will be minimal. 

Analysinq the formation of the shock wave front at the inlet we 
concluded that for the shock wave forminq at the inlet to remain normal to 
the lower wall, inlet should be openinq at the rate: 

v = v , /sin a ( 1 ) 

op sh 7 sk 

where V - velocity of the inlet openinq. 
op 

V , - shock wave velocity in the media, 

sn 

- anqle of the skewed passaqe. 

The velocity is equal to the velocity with which the shock wave surface will 
slide alonq the inlets wall. 

In Fiqs • 2a and 2b results for the modelinq of the Gradual openinq of 
the skewed passaae inlet with the openinq velocity in accordance to equation 
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(2) are shown. As it is obvious from coutour plots of pressure and velocity 



in Figs. 2a and 2b, opening the passage with velocity V calculated from 

op 

eguation (1) assures development of the flow pattern with minimal mixing. 

The shock wave surface in that case remains at all times normal to the walls 
of the passage. 

Conclusions 



Modeling of the gradual opening of the skewed passage revealed the way 
to minimize mixing losses at the passage inlet. The mixing will be minimal 
when the velocity of the opening is matched with shock wave velocity in the 
passage divided by the angle of skewed passage. CXir simulations shows that 
even when conditions for the optimal opening are not satisfied exactly, 
reduction of the rotational mixing could be significant if the opening 
velocity is +15% of the optimal value. 

We can also conclude that a very high opening velocity is reouired in 
order to obtain low mixing loss. The opening velocity of the passage should 
be always higher than the velocity of the shock wave generated in the wave 
rotor passage at the high pressure inlet. Fven substantially skewed passage 
will reouire the opening speed “10% higher than the shock wave velocity in 
the passage which is not easy to achieve for the typical flow conditions in 
the wave rotor. 

Another limitation is that even when this optimal speed of opening is 
achieved at one port, it will not be optimal at other ports; so, 
optimization will not be full. However, mixing losses in the skewed 
passages will be smaller than in rectangular one, if the gradual opening of 
the skewed passage begins from its acute angle. 
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a) Pressure conturs b) Mach number conturs 
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The flow pattern evolution for the passage gradual 
opening with the velocity according to V 0 p £/v 
Time is along the Y - axis (no to scale). 
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Numerical Modeling of Che Non-Steady Thrust Produced by Intermittent Pressure 
Rise in a Diverging Channel. 
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ABSTRACT 



The dynamics of the expansion of detonation 
products in a diverging channel were investigated 
numerically, A two-dimensional unsteady Euler code 
based on the Godunov method was employed* The 
influences of particular features in the expansion 
process , such as the presence of reversed flow and 
propagation of hammer shocks, on the production of 
thrust were examined* Sequentially expansions of 
detonation products were also modeled and it was 
concluded that in order to maintain a high frequency 
periodic mode of operation for propulsion applications 
the channel should be refilled with ambient air after 
each expansion. The influence of the ratio of ambient 
air to detonation product volume on the dynamics of the 
thrust production, and on the impulse generated during 
the expansion, are also reported* 

NOMENCLATURE 

dA x - element of channel internal surface area normal 
to the x-axis 

e - energy per unit volume 

p - pressure 

t - time 

u - component of velocity in x direction 

v - component of velocity in y direction 

x - cartesian coordinate 

y - cartesian coordinate 

e - internal energy per unit mass 

y - ratio of specific heats 

p - density 

subscripts 

0 - conditions at t=0 in combustion products 

ex - external 

L - along inner wall of the channel 
INTRODUCTION 

Early in the development of propulsion engines for 
aircraft, a choice had to be made between engine 
concepts based on nonsteady or on steady gasdynamic 
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processes. Concepts based on steady gasdynamic 
processes were pursued largely because they were 
simpler to analyze and to appreciate conceptually. On 
the basis of the chosen steady concepts, engines were 
built, developed and perfected with time. The 
nonsteady concepts, which looked as promising 
initially, were not developed and remain today in the 
conceptual stage* 

The need for a more thermally efficient small 
engines for various military and civilian 
applications, has led, at intervals, to a 
reexamination of nonsteady concepts ( 1,2,3 ). Today, 
such an examination can be made more thoroughly. For, 
while nonsteady engine concepts were at one time 
extremely difficult to analyze, numerical modeling on 
computers can now provide time-dependent pictures of 
internal flow processes. This can provide efficient 
tools for preliminary design calculations and 
eventually analysis tools for design optimization. 

From the late 50 f s until the early 70 f s, the 
feasibility of an engine operating with intermittent 
detonative combustion was studied at the University of 
Michigan ( 4,5,6 ). Specific impulses over 2100 sec. 
were realized for a single linear shock tube operating 
intermittently with frequencies up to 35 detonations 
per second. Nicholls et al. showed that an engine 
operating intermittently on detonation waves will have 
some advantages over an engine operating on 
deflagration. In general, it will have very high 
specific thrust. The engine will be very simple 
mechanically, and not require precompression of the 
mixture for efficient combustion. 

One of the disadvantages of the engine is that jet 
velocities developed after detonative combustion are 
very high. Efficient operation of the engine for 
propulsion could be realized only at high supersonic 
vehicle velocities (Mach number of about 4). In order 
to reduce the velocity of the out-coming jet, it was 
proposed to use spinning detonation (^6 ) . But the 
spinning detonation process proved to be unsuitable 
for use in an engine because it was unstable. Back et 
al. (jJ) studied the feasibility of using detonative 
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propulsion in Jupiter's atmosphere, and obtained a 
number of results for various types of nozzles. 

Motivated by the potential applications, of 
detonative propulsion, in the present study the 
process of the expansion of detonation products in a 
diverging channel or diffuser of a detonation engine 
were modeled numerically. The model and computer 
program were then applied to calculate repetitive 
firing, and the effect of changing the diffuser 
length. 



THE MATHEMATICAL MODEL AND NUMERICAL SOLUTION 



It is assumed that the jet-propulsion nozzles of 
an engine using detonative combustion can be 
constructed from a number of straight diverging 
channels which are grouped, or clustered together. In 
one cycle of the engine, a small volume of each tube 
is filled with the combustible mixture and undergoes 
detonation. Then the detonation products expand into 
the stationary gas which fills the rest of the tube at 
ambient pressure and temperature. The tube is 
diverging and so the shock wave decays in strength. 
Finally, a weak shock wave leaves the tube, and a 
subsonic flow of gas discharges from the exit of the 
channel. 

It is assumed that the flow is periodic along the 
lines which are a continuation of the walls of the 
channel. This implies that we are modeling the flow 
in a centrally located channel of the cluster. 

It is further assumed that the flow is inviscid. 

The unsteady two-dimensional Euler equations can 
be written in conservation law form as: 



M + 11 + o 

3t 3X 3 Y 
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Here p is the density, u and v are the velocity 
components in the X and Y coordinate directions, p is 
the pressure and y is the ratio of specific heats. 

The energy per unit of volume, e, is defined as 



p( e + 



U 2 + 



-) 



where e = 7 7 - 7 — is the internal energy per unit mass, 

(Y-l)p 

We look for the solution of the system of equations 
represented by Eq (1) in the computational domain as 
shown in Fig 1 in time t, with the following 
conditions at the domain boundaries: 
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a) Solid wall along segment 1-2, 1-3 and 2-4 

The condition on the surface is defined by solving 
the Riemann problem between the point nearest to 
the wall in the domain of integration and its 
mirror image in the direction normal to the wall. 

b) Periodicity between segments 3-5 and 4-6 

By virtue of the central location of the channel 
flow at each grid point along 3-5 is the same as 
at corresponding grid points along 4-6 in Fig 1. 



c) Outflow along segment 5-6 

If the outflow is subsonic at the downstream 
boundary we define only the pressure p out , and 
for all other flow parameters apply the 
continuation condition. That means that u out , 
v ou t anc * Pout are set e< ? ua l to values of u, v 

and p one point ahead of the downstream 
boundary. 

If the outflow is supersonic the continuation 
condition is applied to all parameters at the 
downstream boundary. 

It is assumed that Initially at the time t=0, the 
channel is filled with detonation products from 
segment 1-2 to 1* -2 f . Outside this region of the 
computational domain, the gas is assumed to be at 
ambient conditions . 

The Godunov method has been used to obtain a 
numerical solution of the Eq (1) with the 
described boundary and initial conditions. Details of 
the implementation of the method and boundary 
conditions are given In ( 7) . 

An orthogonal grid is constructed which covers 
the computational domain as shown on Fig 1, It should 
be noted that although the formulation of the problem 
and its solution is two-dimensional, with the boundary 
conditions described above the flow will be 
essentially one dimensional. 

It is assumed that at t*0 , the detonation wave 
has passed through the combustible mixture and the 
products of combustion have uniform properties. Thus 
the time of initiation and propagation of the 
detonation wave through the combustible mixture is 
neglected in comparison with the time for propagation 
through the channel and subsequent expansion of the 
gases . 
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RESULTS AND DISCUSSION 



The history of the flow development in the 
diverging channel is shown in Fig 2. Here the 
velocity and pressure distributions at the center-line 
of the channel are shown at specific moments of time, 
where time t=0 is when the expansion begins. At time 
t=0, approximately 8.5% of the channel volume from the 
left-hand side, is filled with detonation products 
with initial pressure density and velocity having 
values of 

p 0 = 40 atm, p Q ■ 5.2kg/m^, v 0 ■ 0 and u 0 » 0 
respectively. At t=0, the rest of the -channel is 
filled with ambient air with 

p = 1 atm, p = 1.3 kg/m^, v = 0 and u = 0 . 

The length of the channel is 1 meter. These 
conditions will be referred to as Case 1. 

Initially, a shock wave going from left to right 
and a rarefraction wave going from right to left can 
be seen in Fig 2. The shock wave progressively decays 
because the channel is diverging. The rarefraction 
wave propagates rapidly towards the wall at the left 
end of the channel, because the temperatures in this 
region are very high (approximately 2800 K) and the 
channel is converging towards the left end. At a time 
of about 0.15 msec the rarefraction reaches the wall 
at the left end of the channel. From this moment, 
pressure at the left end of the channel rapidly 
decreases from about 40 atm at t 3 0.15 msec to about 
0.3 atm at t = 0.82 msec. Because of this rapid 
depressurization of the region adjacent to the left 
end wall, the pressure on the right hand end of the 
channel becomes higher than that at the left, which 
generates a shock wave going from right to left at 

time t«0.6 msec. This back-going shock wave rapidly 
slows down and then reverses the direction of the flow 
because the channel is converging towards its left end 
wall. 

The negative flow reflects from the left end wall 
at time t«1.4 msec, forming a hammer shock wave which 
moves from left to right and again reverses the flow 
direction. The formation of the hammer shock wave can 
be seen clearly on the velocity and pressure graphs 
for time t = 1.6 msec. Because the shock travels in a 
diverging channel, it weakens rapidly. 

The main shock wave produced by the detonation 
products leaves the channel at the time «0.94 msec, 
and because the channel is diverging in the direction 
of propagation the pressure behind the shock wave 
drops from about 10 atm at time t«*0.1 msec to about 5 
atm, at t = 0.94 msec when it leaves the channel. 
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The mass velocity of the gas drops from about 
1100 m/sec at the beginning of the expansion to about 
500 m/sec when the gas begins to leave the detonation 
tube* The velocity of the gas leaving the tube then 
decreases to -20 m/sec and increases again to about 
50 m/sec when the shock wave reflected from the left 
end wall reaches the exit of the channel. This 
concludes the significant events which occur in the 
diverging channel as a result of the expansion of the 
detonation products. At later times, rarefraction and 
pressure waves continue to form and travel in the 
channel, but they are much weaker and their influence 
on the thrust produced is very small. 

In Fig 3 thrust as function of time is plotted 
for the case described in the graphs given in Fig 2. 
Thrust was calculated for the case of a stationary 
tube with external pressure equal to the ambient 
pressure using 

F - / L (p L - P ex )dA x (2) 

It can be seen in Fig 3 that most of the thrust 
is produced in the 2 msec following the beginning of 
the expansion of the detonation products. Integration 
of thrust in time shows that 90% of the impulse is 
produced in the first 2 msec and 10% in the following 
7 msec* 

The peculiarities of the thrust curve can be 
understood in relation to the wave pattern in the 
tube. The first change in slope at time t«0.8 msec 
corresponds to the time when the front of the main 
shock wave reaches the end of the tube. Because the 
pressure behind the shock front drops rapidly, the 
thrust decay increases when the main shock front 
leaves the channel. At time t«1.5 msec the thrust 
increases as a result of the hammer shock wave 
reflecting from the left wall. 

When the main process of expansion ends, at 
time t«6 msec, the pressure in the tube is about 1.1 
atm and the average density is about 0.43 kg/m^. This 
corresponds to an average temperature in the channel 
of about 935 K. Since the pressure in the tube is 
close to ambient, the heat exchange will be by natural 
convection, and it will take a relatively long time 
for the gas remaining in the tube after the expansion 
to cool. 

In order to examine what will be the wave pattern 
and thrust when a second detonation wave expands in 
the channel immediately following the expansion of the 
first detonation wave (6 msec after the beginning of 
the process of expansion of the first detonation 
wave), conditions were simulated numerically in the 
following way. Results of the calculation for the 
first detonation wave expansion were used as the 
initial distributions of flow parameters. Then, in 
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the 8.5% of the channel volume adjacent the left end 
wall the pressure, density and velocity were given the 
values 

p Q =* 40 atm, p Q = 5,2 kg/ m 3 and v 0 = u Q * 0 
respectively at time t=0. These conditions will be 
referred to as Case 2. 

In Fig 4, the thrust vs time is shown for the 
second detonation wave expansion* Since the 
detonation products expand in a low density hot gas, 
the expansion goes much faster than in the previous 
case. At time t«0.8 msec from the beginning of the 
second detonation wave expansion, the thrust drops to 
zero. This can be compared with a thrust of 
approximately 6*10^N/m at the corresponding time shown 
in Fig 3, for expansion into ambient air. Then, 
because the large volume of the channel has static 
pressure lower than ambient, thrust becomes negative 
and, as in the previous case, a back-going shock wave 
develops. In Fig 5, the development in time of this 
recursive shock can be followed. It is observed that 
the negative flow velocities which develop are twice 
as large as those for the preceding expansion into 
ambient air. And at time t*2 msec the flow throughout 
the tube is reversed. Reflection of the 
left-traveling shock wave from the left end wall at 
time t»2 msec the flow throughout the tube is 
reversed. Reflection of the left-traveling shock 
wave from the left end wall at time t«2 msec produces 
a very strong hammer shock which reverses the thrust 
and the flow direction at time t«2.1 msec. Because of 
the large negative thrust, the total impulse which is 
produced in Case 2 is about 15% smaller than the 
impulse of the detonation products expanding in 
ambient air. 

Another notable effect of expansion into the 
products of the previous cycle, is the substantial 
rise of temperature which occurs in the channel after 
the second expansion. When the process of the second 
expansion ends the average temperature in the channel 
is approximately 1550 K. This is 615 K higher than 
the final temperature after expansion into ambient 
air. 

In order to examine how thrust and total impulse 
are influenced by the volume of air contained at 
ambient conditions in the channel before the 
expansion, two additional test cases were modeled. In 
Case 3, the volume of ambient air in the channel to 
the right side of the interface with the detonation 
products is taken to be approximately 2.5 times larger 
than the volume of the detonation products. This 
requires a channel about 50% shorter than the original 
channel, for which the ratio of the volumes of ambient 
air to detonation products was about 11. In Case 4, 
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there was no air in the channel when expansion began, 
which implies that the channel in this case has 
only about 25% of the original length* Therefore, 
between Cases 3 and 4 and Case 1 only the length of 
the channel was changed* The geometry of the channel, 
initial parameters of the detonation products were the 
same in Cases 1, 3 and 4. 

In Fig 6 the thrust is shown as a function of 
time for Cases 1, 3 and 4. It can be seen that Case 1 
(solid line) where the channel is fully extended, 
gives consistently higher thrust and that the thrust 
decreases when the volume of ambient air in the 
channel (or channel length) decreases; Case 3 (chain- 
dotted line) and Case 4 (dashed line) have regions 
where thrust is negative* Comparison of the total 
impulse generated shows that the extension of the 
channel is very beneficial for producing impulse. The 

impulse produced in Case 4, where the tube is filled 
with detonation products is 0 . 162 • lCpN/M *sec • When 
the channel is extended so that only 28.6% of the 
channel volume is initially filled with detonation 
products (Case 3), the impulse increases 28%. An 
additional extension of the channel so that only about 
8.5% of the tube volume is initially filled with 
detonation products (Case 1), leads to total impulse 
of 0.2784*10^ N/M*sec, which is a 72% increase 
over Case 4. 

CONCLUSIONS 

Numerical modeling of the process of expansion of 
detonation products in a diverging channel revealed a 
flow pattern rapidly changing in time with multiple 
shock and rarefraction waves. One of the most 
interesting features was the occurrence of a 
backwards-traveling shock wave, which developed as a 
result of over-expansion in the channel. The shock 
wave first reversed the flow direction (and in some 
cases developed a negative thrust) and then, 
reflecting from the left end wall as a hammer shock, 
reversed the flow direction again, increasing the 
thrust . 

Modeling of the second detonation wave expanding 
into the products of the previous expansion, showed 
that this arrangement leads to: a) significant 

negative thrust; b) increases of the gas temperature 
after expansion; c) losses of 15% of the total 
impulse. This leads to the conclusion that the 
channel should be filled with ambient air after each 
detonation and expansion. Since the process of 
expansion takes only about 0.006 sec, the 
time-averaged thrust which can be generated by a 
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single channel will be limited by the rate at which 
the channel can be re-charged with air and mixture, 
and refired. For very low initial velocities in the 
channel (consistent with the initial conditions 
specified in the modelling), the frequency of firing 
is limited to about 20 detonations per second. This 
corresponds to 0.044 sec for recharging at gas 
velocities of about 23 m/sec. Investigation of the 
influence of the channel length on the impulse and 
thrust production showed that increases in the volume 
of ambient air in the channel ahead of the detonation 
products, which was obtained with an increase in 
channel length, led to significant increases in the 
thrust and impulse generated during the expansion. 

The total impulse generated by an expansion in the 
channel in which 8.5% of the total volume was filled 
initially by the detonation products, was found to be 
72% greater than in the case of a short channel filled 
only with the same volume of detonation products at 
the same initial conditions. 

The diverging channel geometry provided the 
simplest geometry with area change to model using the 
unsteady 2D Godunov-Euler code. Variations in the 
geometry can now be made in order to examine the 
unsteady thrust performance. 
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Figure 1. The Computational Domain. 
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Figure 2. 



The Evolution in Time of the Detonation 
Products Expansion. 
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Figure 3. The Evolution in Time of the Thrust 

Developed by the Expansion of the Detonation 
Products . 
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Figure 4, The Thrust vs. Time Curve for the Second 
Sequential Expansion of the Detonation 
Products . 
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Figure 5. Formation and Evolution of the Negative 

Flow Velocities for the Second Sequential 
Expansion of the Detonation Products. 
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Figure 6. The Thrust vs. Time Curves for Expansion 
in Channels of Different Lengths. 

Solid Line - Longest Channel 
Dashed Line - Shortest Channel 



APPENDIX C 



Numerical Techniques for Wave Rotor Cycle Analyses. 
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INTRODUCTION 

Wave Rotors are devices which use unsteady 
pressure waves to effect direct energy exchange between 
two gases. They offer the potential for significant 
Improvements In air-breathing engine propulsion cycles 
through their capability of withstanding higher 
temperatures and pressures than- present-day turbines. 
Other diverse applications are recounted In Ref (1). 

In Its most basic form, a wave rotor Is a drum 
with axial or helical (staggered) passages arranged 
around the periphery. This single drum-like rotor 
replaces separate compressor and turbine components for 
gas turbine applications. The compression and 
expansion of the two fluids occurs In the passages as a 
result of shock tube like processes. In a typical 
configuration, combustion products at high temperature 
and pressure give up energy to the ’driven' fluid 
(usually air) through the action of time-unsteady 
compresslon/shock waves. The combustion products, (hot 
'driver' gases), are In turn expanded and exhausted 
from the rotor, the available work of expansion being 
utilized to Induce a fresh charge of air onto the 
rotor. Careful sequencing of the passage ends past 
stationary Inlet and outlet ports In valve plates on 
either side of the rotor creates a cyclic Internal wave 
pattern In the wave rotor component. The high 
temperature capability of the device is a direct 
consequence of the repetitive processing of both cold 
and hot fluid alternately in the same rotor. 

Estimation of the aero- thermodynamic performance of 
these devices hinges on calculations of the unsteady 
energy exchange in the wave rotor component. Numerical 
simulation of the unsteady wave processes can generally 
be carried out on a one-dimensional basis. However, 
the complex pattern of flow discontinuities and wave 
interactions known to exist in wave rotors call for 
numerical methods which solve the nonlinear, hyperbolic 
system of governing equations without relying on either 
artificial viscosity or special treatment of 
discontinuities. Described in this paper are two such 
techniques, the Godunov method and a Random Choice 
method (Glimm's method or Piecewise Sampling method). 



NUMERICAL FORMULATION 

An essential component of both the Random Choice 
method and the Godunov method Is the solution of a 
sequence of local Riemann problems. A Rlemann problem 
Is defined by setting up an Initial value problem for 
the equations of motion in Eulerian form. The 
unsteady, two-dimensional Euler equations can be 
written in vectorized, conservation law form as: 

U t + (F(U)] X + [G(U) Jy - 0, where U ■ |p pi pv e] T , 

P(U) - [ pu (p + pu 2 ) puv (e + p)u]T, and 

G(U) -«■ [ pv puv (p + pv 2 ) (e + p)v] T . 

Here, p is the density, u and v are the velocity com- 
ponents in the x and y coordinate directions, and p is 
the pressure. Subscripts indicate partial differenti- 
ation with respect to that independent variable. The 
energy per unit volume, e, is defined by: 

e - p( c + (u 2 + v 2 )/2), where t - p/(y -1 ) p 

is the Internal energy per unit mass. 

The Riemann problem is intrinsically 
one-dimensional in nature and is defined accordingly. 
For the one space variable case, thus, U - U(X,t), and 
the initial value problem is set up by specifying 
initial conditions which consist of intervals over 
which data are constant, separated by jump 
discontinuities. If the time step is chosen to be 
sufficiently small, the waves propagating with finite 
speeds from adjacent discontinuities remain within 
their respective spatial cells and do not intersect. 
This sequence of solutions from adjacent Riemann 
problems is then pieced together to obtain the whole 
solution at each succeeding time step. 

The Random Choice method (RCM) for gasdynamics is 
the outcome of a constructive existence proof of 
solutions to systems of nonlinear hyperbolic 
conservation laws presented by Gliram Ref (2), and its 
development into an effective numerical tool by Chorin 
Refs (3,4). The solution is advanced in time by a 
method that Includes a solution of Riemann problems as 
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described above and a sampling procedure, which 
chooses values from a representative point of the 
exact solution to the local Rlemann problem. The 
sampling procedure is due to van der Cor put Ref (5) 
and generates equidistributed sequences Ref (6). The 
technique eliminates numerical diffusion and flow 
discontinuities are computed with infinite resolution, 
i.e. , there is no smearing. 

The Godunov method is a finite difference method 
which computes the Rlemann problem solutions as a 
first step, and using the fluxes thus obtained, 
advances in time by solving the first-order finite 
difference approximation of the Euler equations. 

Coding details for the RCM, and the Godunov 
method are given in Refs (4,7) and Refs (7,8) 
respectively. In the following section, examples of 
the applications of these techniques to wave rotor 
devices are given. 

EXAMPLE APPLICATIONS 

Fig 1 shows a wave diagram (a one-dimensional 
space-time plot) for a wave rotor device configured to 
operate as a turbine alone. Air at a pressure higher 
than ambient enters the rotor through the inlet port, 
generating a shock wave with a slower moving contact 
discontinuity behind it at point 'a'. The shock 
reflects off a wall boundary at point f b f , crosses the 
Incoming interface at point r c' (bringing it to a 
halt), and reaches the inlet side again at point 'd', 
whereupon the inlet port is dosed. The gases in the 
rotor passages are now in an essentially quiescent, 
high pressure state. At point r e ? , the outlet port is 
opened to a lower pressure, a rarefaction fan being 
generated in the process. The expansion waves travel 
to the left, reflect off the solid boundary and 
propagate back to the outlet side. The port is dosed 
when conditions in the passages are essentially the 
same as at the beginning of the cycle. 




Fig. 1 Wave Diagram Computed by 1-D Random Choice 
Method. 

S— Shock; RS— Reflected Shock; R— Rare faction 
Fan; RR-Reflected Rarefaction; I-Interface 



The entire wave diagram was generated by the 
one -dimensional sampling method with properly 
Implemented boundary conditions. Fig 2 shows a 
sequence of the propagation of the shock wave and 
contact discontinuity in time. The perfect resolutic 
of the discontinuities is noteworthy, even when a vei 
small density change occurs across the interface, as 
in this case. 




Fig. 2 Sequence Showing Shock and Interface Movemeni 
S - Shock; I - Interface; RS - Reflected 
Shock; N - Times tep 

The transient process of the rotor passage end 
opening or closing creates losses which affect the 
performance of the device. The gradual (as opposed 1 
instantaneous) opening/closing of the passages is 
essentially a two-dimensional flow phenomenon and is 
modelled using the 2-D Godunov code. Fig 3 shows a 
sequence of pressure contours for the gradual opening 
case which clearly shows the significant effect the 
dynamics of the passage opening has on the flow 
pattern in wave rotor devices. The actual shock 
formation for the modelled case is seen to occur 
approximately halfway into the passage as opposed to 
the Instantaneous formation generally assumed for the 
ideal case. 

DISCUSSION 

The two techniques described here are powerful 
tools for the design of wave rotor devices. The 
example application of the one -dimensional Random 
Choice method deals with a fairly elementary wave 
diagram, but serves to illustrate its potential and 
suitability for designing complex gas turbine engine 
cycle type wave diagrams. The Godunov method is mor€ 
amenable to extension to multidimensional analysis 
required for loss calculations, while fulfilling the 
requirements outlined in the introduction. 
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Fig. 3 Contour Plots Generated by 2-D Godunov Method 
for Gradual Opening of Wave Rotor Passage End. 
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APPENDIX D 



Local Cell Orientation Method - Illustration of the Method for the Solution 
with an Oblique Shock. 
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Local Cell Orientation Method. 



Illustration of the method for the solution with oblique shock 



INTRODUCTION 



The Finite Volume Schemes for numerical solution of the Euler equa- 
tions can be split into two main steps: 

a) Calculation of the fluxes on the cell edges; 

b) Updating the values at the center of the cell , using the 

calculated fluxes, to satisfy the Euler conservation laws. 

Usually the accuracy of the flux calculation determine the accuracy 
of the numerical modeling Ref.1. In the region of the shock wave the 
accuracy of the flux calculation is significantly reduced and in some 
methods special care should be taken in order to reduce the pre-shock 
and after-shock oscillations .When the upwind methods are used the 
fluxes are usually calculated using the one dimensional wave propagation 
problem (i.e. Riemann problem). In this case the assumption of the 
one-dimensionality lead’s to large computational errors and significant 
shock smearing in the region of the oblique shock. 

In the current study it will be demonstrated that the accuracy of 
the numerical modeling could be significantly improved by proper Local 
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Cell Orientation (SELCO Method). We applyed the SELCO Method only 
with the Godunov scheme for the solution of the Euler equations, but 
the idea of the method is general enough to be used with other finite 
volume methods. 



THE SELCO METHOD. 



For the Godunov method the numerical fluxes are calculated using 
the solution of the Riemann problem. First the solution is assumed to be 
piece-wise constant in every cell. Then, to calculate the flux for every 
cell edge the Riemann problem is solved between the cells adjacent to 
the edge. Solution of the Riemann problem is obtained in the direction 
normal to the cell's edge. The equations are then integrated in the 
each cell for one time step using the finite difference approximation of 
the Euler equations and the fluxes at the cell's edges. A more complete 
description of the Godunov method can be found in Ref.1. 

In the two dimensional Godunov method it is assumed that the fluxes 
could be determined using the solution of the one dimensional Riemann 
problem if the CFL number used is smaller than 0.5. In the following 
example it will be demonstrated that that assumption lead's to consider- 
able errors for the solutions with oblique shock waves. 

As a test case the supersonic flow over the 22.3° wedge was simu- 
lated. The grid and the computational domain is shown on the Figure 1. 
The system of the two dimensional Euler equations as in Ref.1, was 
solved by time marshing in the computational domain shown in Fig.1 for 
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time t * oo with the following conditions at the domain boundaries: 

a) Inflow of the air at M=2.2 along segment 1-2; 

b) Outflow along segment 3-4; 

(the flow along this segment will be supersonic and 
the continuation condition is appli-ed to all parameters 
at the downstream boundary) 

c) Solid wall along segments 1-3 and 2-4. 

In the Figure 2 (a,b,c,d) results of the flow simulation using the 
Godunov method in the computational domain shown on Fig.1 (squares) 
are compared with the analytical solution for the supersonic flow over 
the wedge (solid lines) Ref. 2. Comparison is done for the data on the 
lower surface of the channel shown on the Fig.1. It can be concluded 
from this comparison that only pressure coefficient is calculated accu- 
racy by the Godunov method. Especially notable the error in entropy. 
Exact entropy change on the shock wave is 2 times smaller than that 
predicted by the Godunov method. The isomach lines are shown for the 
same simulated case, in Figure 2 (e) . The shock is smeared over the 
large area of the mesh (up to 5-6 points), but the shock's angle is 
simulated correctly. 

In order to study the source of errors for this type of problems the 
supersonic flow over the wedges with different angles was simulated. It 
was found that the pressure was predicted accuratly for the all simu- 
lated cases. At the same time the error in entropy increases when the 
shock wave obliquity increases. Dependence of the accuracy of the 
shock simulation on the shock obliquity is not surprising since the basic 



53 



PAGE 4 



assumption about independance of the flux calculation for the intersect- 
ing edges of the cell, for t < 0.5*CFL, became incorrect in the vicinity 
of the shock. Spliting the flux calculation for each cell on four one 
dimensional Riemann problems for each of the cell's edges would not 
lead to the large errors in the regions of monotonic flow. It will also 
give accurate approximation in the region of the oblique shock if the 
shock wave is parallel to one of the cell edges. The last statement will 
be illustrated by the following numerical example. 

Let's solve the problem of the supersonic flow with M =2.2 over the 
22.3° wedge, as in previous example, on the grid shown in Figure 3. 

The grid in Fig.1 is skewed uniformly so that vertical lines of the grid 

are forming a 51.7° angle with the positive direction of the x-axis. An 
angle of 51.7° corresponds to the shock wave angle calculated analyt- 
ically Ref. 2. Because the grid is skewed the oblique shock waves will 

be parallel to the vertical grid lines of the mesh shown in Fig. 3. 
Results of this test case, in the form of the distribution of pressure 
coefficient, entropy, density and velocity, are compared to the analyt- 
ically calculated values (solid lines) in Figure 4 (a,b,c,d). In Fig. 4 
(e), the isomach lines for this simulation are shown. It can be con- 
cluded from the results presented in Fig. 4, that the flow is modeled is 
this case with very high accuracy. Shock wave is resolved on one grid 
cell and entropy jump is calculated exactly. All that is result of the 
proper cell orientation towards the shock wave. 

Skewing the entire grid can help to accuratly resolve only one 
shock wave, and only on the condition that the shock wave angle is 
known prior to the solution, so it could not be applied effectivly. But 
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if skewing of the cell can be done locally in the region of the oblique 
shock it can improve the accuracy of the shock wave simulation. 

The proposed SELCO method does the cell reorientation locally. The 
method consist from the following steps: 

1 ) I ntegration of the Euler equations by the-Godunov method; 

2) Define the approximate shock location and the shock angle using 

the expressions for oblique shocks Ref. 2. 

(Based on the fact that pressure is calculated accuratly 
by the Godunov method); 

3) Rotate the cell edges that are directed along the shock about 
their middle points on the shock wave angle. So, after rotation, 
two of the cell edges will be parallel locally to the shock wave 
surface (see Figure 5); 

4) Calculate the fluxes on the new edges of the cell; 

5) Integrate the Euler equations in the new cell. 

All these additional steps do not add much of computational work, 
because the cell reorientation should be done only in the vicinity of the 
shock wave. If the shock could be resolved on one grid point, only one 
cell would be transformed. Our experience has shown that there is no 
need to extend the SELCO procedure on the cells ahead and behind the 
shock wave surface. 

In Figure 5 (a,b,c,d) results are shown for the same flow condition 
as in previous cases: supersonic flow with M =2.2 over the wedge of 
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22.3°. Calculations were done on the grid shown in Figure 1 and in 
this case the SELCO method was applyed. It can be concluded from the 
results presented in Fig. 5 that the accuracy of the shock wave model- 
ing using the SELCO method approach that demonstrated in Figure 4 for 
the completely skewed grid, and is superior to the accuracy of the 
standard Godunov method. The isomach lines for the case where the 
SELCO method was applied are presented in Figure 5 (e). The shock 
tickness in this graph is minimal and much more improved compared to 
the simulation by the original Godunov method shown in Figure 2 (e) . 



CONCLUSIONS. 



It has been demonstrated that inaccurate modeling of the oblique 
shock waves produced by the Godunov method is the result of the 
obliqueness of the shock wave with respect to the edges of the cells of 
the computational grid covering domain of integration. It is also shown 
that only when the shock surface is parallel to the two opposite edges 
of the cell the oblique shock can be accuratly calculated. 

A new method of the Local Cell Orientation ( SELCO Method) 
is proposed in order to allow local reorientation of the cells in the 
vicinity of the shock waves. The efficiency of the SELCO method is 
demonstrated for the simulation of the oblique shock waves in the 
supersonic flow using Euler equations. 

Although the new SELCO method was demonstrated with the Godunov 
scheme it will be effective in applications to the other upwind methods 
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that use the finite volume formulation. 
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FIGURE CAPTIONS. 

1. The Computational Domain 

2. Computation by the Godunov method . =2 . 2. 

Wedge Angle is 22.3°. Continuous lines 
show the exact solution for - a,b,c and d. 

a) Surface Pressure Coefficient. 

b) Surface Entropy. 

c) Surface Density. 

d) Surface Mach Number. 

e) Isomach Lines. 

3. The Computational Domain with the Skewed Grid 

4. Computation by the Godunov Method on the Skewed Grid Shown 

in Figure 3. M_ oo =2.2. Wedge Angle is 22.3°. 

Continuous lines show the exact solution for - a,b,c and d. 

a) Surface Pressure Coefficient. 

b) Surface Entropy. 

c) Surface Density. 

d) Surface Mach Number. 

e) Isomach Lines. 

5. The Local Cell Reorientation by the SELCO Method. 

6. Computation by the SELCO Method . M_„=2 . 2 . 

Wedge Angle is 22.3°. Continuous lines 
show the exact solution for - a,b,c and d. 

a) Surface Pressure Coefficient. 

b) Surface Entropy. 

c) Surface Density. 

d) Surface Mach Number. 

e) Isomach Lines. 
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The Computational Domain. 
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2. Computation by the Godunov method . M_ oo =2. 2. 
Wedge Angle is 22.3°. Continuous lines 
show the exact solution for - a,b,c and d. 
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2. Computation by the Godunov method.M co =2.2. 
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a) Surface Pressure Coefficient. c ) Surface Density. 

4. Computation by the Godunov Method on the Skewed Grid Shown 

in Figure 3. M oo =2.2. Wedge Angle is 22.3°. 

Continuous lines show the exact solution for - a,b,c and d. 
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The Local Cell Reorientation by the SELCO Method. 
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a) Surface Pressure Coefficient. c ) Surface Density. 

6. Computation by the SELCO Method . M^=2. 2. 

Wedge Angle is 22.3®. Continuous lines 
show the exact solution for - a,b,c and d. 
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